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Introduction 


This paper describes the implementation of optimization techniques based on control theory for airfoil design. 
In our previous work in the area [6, 7, 19, 11] it was shown that control theory could be employed to devise 
effective optimization procedures for two-dimensional profiles by using either the potential flow or the Euler 
equations with either a conformal mapping or a general coordinate system. We have also explored three- 
dimensional extensions of these formulation* recently [10, 9, 20]. The goal of our present work is to demonstrate 
the versatility of the control theory approach by designing airfoils using both Hicks- Henne functions and B- 
spline control points as design variables. The research also demonstrates that the parameterization of the 
design space is an open question in aerodynamic design. 

Formulation of the design 
problem as a control problem 


in this research control theory serves to provide computationally inexpensive gradient information to a standard 
numerical optimization method. For flow about an airfoil or wing, the aerodynamic properties which define 
the cost function are functions of the flow-field variables (w) and the physical location of the boundary, which 
may be represented by the function T , say. Then 

/ = liw,T) 


and a change in T results in a change 

U 

ow OS 

in the cost function. Each term in (1), except for 6w, can be easily obtained. Finite difference methods 
evaluate the gradient by making a small change in each design variable separately, and then recalculate both 
the grid and flow-field variables. This requires a number of additional flow calculations equal to the number of 
design variables. Using control theory, the governing equations of the flowfield are introduced as a constraint 
in such a way that the final evaluation of the gradient does not require multiple flow solutions. In order to 
achieve this, Sw must be eliminated from yl). The governing equation R and its first variation express the 
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dependence of u and T within the flowfield domain I) 
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Next , introducing a Lagrange multiplier v. we have 
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Choosing i to satisfy the adjoint equation 
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the first term is eliminated, and we find that the desired gradient is given by 
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The advantage is that (4) is independent of 6u\ with the result that the grad.ent of / with respect to an 
arbitrary number of design variables can be determined without the need for additional flow-field evaluations 
The main cost is in solving the adjoint equation (3). In general, the adjoint problem is about as complex as 
a flow solution. If the number of design variables is large, the cost differential between one adjoint solution 
and the large number of flowfield evaluations required to dete-mine the gradient by finite differencing becomes 
compelling. Once equation (4) is obtained. Q can be fed into any numerical optimization algorithm to obtain 
an improved design. 


Issues of importance for design problems 

The develop nent of aerodynamic design proceedures that employ an adjoint equation formulation is currently 
being investigated by many researchers. These methods promise to allow computational fluid dynamics meth- 
ods to become true aerodynamic design methods. References [1. 2. 12, 13, 16 15. 14. 21. 17, 22] represent a 
partial list of recent works that reflect this developing field. However, as is the case in any new research field, 
many questions remain Probably the most salient issues of concern are the following: 

1. Discrete vs. continuous sensitivities 

2. Choice of optimization procedures 

3. Treatment of geometric and aerodynamic constraints 

4. The level of coupling between design and analysis 

5. The parameterization of the design space 

These topics, which were outlined in a proposal to the P S. Airforce in 1991 [8], still require intensive inves- 
tigation in the up coming years. With regards to the first item, it is historically interesting that Jameson in 
! n 88 [6] first developed the equations for a continuous approach to transonic airfoil and wing design, whereby 
the steps represented by equations (1-4) were applied to the governing differential equations. The resulMng 
adjoint differential equations with the appropriate boundary conditions may then be discretized and solved 
in a manner similar to thf.t used for the flow solver. One may alternatively derive a set of discrete adjoint 
equations directly from th •> discrete approximation to the flow equations by following the procedure outlined 
in equations (1-4). The resulting discrete adjoint equations are one of the possible discretizations of the con- 
tinuous adjoint equations. T his alternative is mentioned in [7], but was not adopted in that work because 
of the complexity of the resulting discrete adjoint system. It has the advantage that the resulting discrete 
gradients are the true gradients of the discrete model of the flow equations, and it has found favor in the work 
of Taylor et al. [16. 13, 12] and others [1. 2]. It seems that both alternatives have some advantages. The 
continuous approach gives the researcher some hope for an intuitive understanding of the adjoint system and 
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its related boundary conditions. It also permits an easy recycling of the solution methodology used for the 
flow solver. The discrete approach, in theory, maintains perfect algebraic consistency at the discrete level If 
properly implemented it will give gradients which closely match those obtained through finite differences. The 
continuous formulation produces slightly inaccurate gradients due to differences in the discretization. H owe vp r, 
these inaccuracies must vanish as the mesh width is reduced. A drawback of the discrete approach is that, 
once the large matrix problem of equation (3) is defined it becomes more difficult to recycle the flow solution 
algorithm on the resulting linear algebra problem. It is subject, moreover, to the difficulty that the discrete 
flow equations often contain nonlinear flux limiting functions which are not differentiable, it also limits the 
flexibility to use adaptive discretization techniques with order and mesh refinement, such as the h-p method, 
because the adjoint discretization is fixed by the flow discretisation 
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Design variables 


It turns out that the determination of items (2-4) in the above list strongly hinge on the choice for (5). In 
Jameson's first works in the area [6, 7. 10], every surface mesh point was used as a design variable in three- 
dimensional wing design cases this led to as many as 4224 design variables [10] The use of the adjoint met hod 
eliminated the unacceptable costs that such a large number of design variables would incur for traditional 
finite difference methods. If the approach were extended to treat complete aircraft configurations, at least tens 
of thousands of design variables would be necessary. Such large numbers of design variables preclude the use 
of descent algorithms such as Newton or quasi-Newton approaches simply because of the high cost of matrix 
operations for such methods. The limitation of using a simple descent procedure such as steepest descent 
has the consequence that significant errors can be tolerated initially in the gradient evaluation. Therefore 
such methods favor tighter coupling of the flow solver, the adjoint soiver, and the overall design problem to 
accelerate convergence. Ta’asan et al [21, 14] took advantage of this by formulating the sign problem as a 
one shot procedure where all three systems are advanced simultaneously. Choosing such * jesign space suffers 
from admitting poorly conditioned design problems. This is best exemplified by the case where only one point 
on the surface of an airfoil is moved, resulting in a highly nonlinear design response. In his original work 
Jameson [7, 9] treated this poor conditioning of the design problem by smoothing the control (surface shape) 
and thus removing the problematic high frequency content from the advancing solution shape. 

Further, arbitrarily freeing all surface points as part of the design makes it difficult to enforce geomet- 
ric constraints, such as maintaining fuel volume, without resorting to constrained optimization algorithms 
Nevertheless, this choice for variables remains attractive since it admits the greatest possible design spa'-e. 

Ray Hicks and others [5,4, 18] have in the past parameterized the design space using sets of smooth functions 
that either generate or perturb the initial geometry. Bv using such a parameterization it is possible to work 
with considerably fewer design variables than the choice of every mesh point. Hicks initially adopted the 
approach because his work in the past exclusively used finite difference based gradient design methods that 
are inherently intolerant of more than a few dozen design variables for large problems. However, since these 
early methods had to content themselves with at most ?. hundred design variables, great freedom in the choice 
of the design algorithm was possible. One simple choice of design variables for airfoils suggested by Hicks and 
Henne [4], have the following ”sine bump” form: 



Here l\ locates the maximum of the bump in the range 0 < r < 1 at x = fi, since the maximum occurs when 
x Q = where a = log log 1 1 i, or a\ogt\ ~ !og^. Parameter 1 2 controls the width of the bump. 

When distributed over the entire chord on both upper and lower surfaces, these analytic perturbation 
functions admit a large possible design space. They can be chosen such that symmetry, thickness, or volume 
can be explicitly constrained, thus avoiding the use of expensive constrained optimization algorithms to address 
geometric constraints. Further, particular choices of these variables will concentrate the design effort in regions 
where refinement is needed, while leaving the rest of the airfoil section virtually undisturbed. The disadvantage 
of these functions is that they are not orthogonal, and there is no simple way to form a basis from these functions 
which is complete for the space of continuous functions which vanish at x = 0 and x — 1 . Thus, they do not 
guarantee that a solution, for example, of the inverse problem for a realizable target pressure distribution will 
necessarily be attained. Nevertheless, they have proved to be quite effective in realizing design improvements 
with a limited number of design variables. The design process that usee these basis spaces can be accelerated 
toward convergence by either tighter coupling of the individual design elements, as was the case when using 
the mesh points themselves, or through the use of higher order optimization algorithms Finally they have 
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one last advantage over using the mesh points; there is no need to smooth the resulting solutions as the design 
proceeds when using Hicks- Henne functions, sine*' by construction higher frequecies are not admitted and thus 
the design spaces are naturally well posed 

Another set of design variables that has recently achieved favor is the us*- of B-spline control points In this 
approach the geometry is defined by B-spline curves or surfaces while their control points become the design 
variables. Like the Hirks-Henne functions. B-splines allow for a greatly reduced number of design variables, 
and thus permit the use of. say. a quasi-Newton design procedure. If the upper and lower surfaces of an airfoil 
are separated, the method easily admits camber or thickness constraints explicitly within the design spao- 
Further, local control is also possible by choosing only a limited number of control points as active design 
variables. This method in practice seem- to have an advantage over the Hicks Henne functions in that a more 
complete basis space of admitted airfoils is permitted with the same number of design variables Finally since 
these curves and surfaces are the natural entities u.^d in (’AD envirornents, they provide a straightforward 
way of integrating (’AD and aerodynamic design 


Numerical Experiments 

In this short paper, a study of both Hicks- Henne and B-spline control points is explored to begin to address 
the open issue presented in item (5) above. For the details of the development and implementation of the 
full design procedure, the reader is referred to our previous works listed earlier. Here the two-dimensional 
airfoil design problem is accomplished by combining an Euler flow solver, a corresponding adjoint solver, and 
an unconstrained quasi-Newton optimization algorithm that uses BFGS updates [3, 11]. lh*° method has 
beenj demonstrated to be successful in the design of various airfoils at transonic conditions and with various 
objective functions. 

Figure (1) shows a design case where 25 Hicks-Henne burnps are used as design variables for both the upper 
and lower surfaces. The bumps are arranged with high clustering toward the leading and trailing edges to give 
adequate design fidelity in these sensitive areas. Figure (la) shows the initial conditions for an inverse design 
problem where a target pressure distribution for an RAE 2822 airfoil at a = 1°. M = 75, and C\ ~ .6030 is 
specified and a NACA 0012 airfoil is prescribed as a starting airfoil. Both the target and the initial airfoil show 
strong shocks in their respective solutions. After 5 design iterations significant progiess is seen in matching the 
target pressures. Moreover, the cost function, which is measured as an integral of the square of the differences 
between the two pressure distributions, has dropped from 68.556 to 4.247. The proper location for the shock 
has already been achieved. By 10 iterations the cost function has been further reduced to 0.702. and after 30 
iterations no differences are visible in the solution where the cost function is only .010. 

Figure (2) shows the same design problem treated with two B-spline curves for the upper and lower surfaces. 
Each curve was defined by !8 control points with a degree of 4 and was chosen to represent the airfoil shape 
directly. The y ordinate of each control point was used as a design variable. A close approximation to a NACA 
0012 airfoil is again used as the initial point with the same target pressure distribution of an RAE 2822 airfoil 
at ft = 1°, A/ = .75, and Ci = .6030 being specif 1. It is seen that the final cost function of .005 is lower than 
that of the final solution for the Hicks-Henne functions. However, not only did this solution take 50 design 
iterations as opposed to 30 for the first example, but the interim solutions displayed in Figures (2b and 2c) 
show gross oscillations in the design. The conclusion that may be drawn from this is that the use of B-spline 
control points localizes the design variables to a point that they act more like using the mesh points themselves 
as design variables. Thus the high frequency information contained in the control can easily cause oscillations 
in the design, with the result that there is a greater degree of nonlinearity and a poorer conditioning of the 
design space as compared to that for the Hicks-Henne functions. A prudent step would be to use a low pass 
filter when using B-spline control points as design variables to smooth the control in a manner similar to that 
employed by Jarneson [6, 7. 10]. Attempts to increase or decrease the degree and the number of control points 
for the B-spline curves did not prove to have favorable effects. 
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la Initial point 
Cost function = 68 556 
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lb Design after 5 iterations 
Cost function = 4 247 



Figure 1: Inverse design from a NACA 0012 to an RAE 2822: Hicks-Henne design variables 
Af - 0.75. Ci = 0.6031. C d = 0.0046, o = 1.000. 50 Design Variables 
— . x Initial Airfoil NACA 0012 
- - + Target C p RAE, M = 0.75 

193x33 Mesh, Euler Solution 
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2a: Initial point 
Cost function = 68.605 



2b: Design after 5 iterations 
Cost function = 16.500 



Figure 2: Inverse design from a NACA 0012 to an RAE 2822 E-spline control point design variables 
M = 0 75. C { = 0.6030. C d = 0.0045. a = 1. 000. 50 Design Variables 
— x Initial Airfoil: NACA 0012 

. + Target C v RAE. M = 0.75. 

193x33 Mesh, Euler Solution 
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